Multi-critical point in a diluted bilayer Heisenberg quantum antiferromagnet 
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The S — 1/2 Heisenberg bilayer antiferromagnet with randomly removed inter- layer dimers is 
studied using quantum Monte Carlo simulations. A zero-temperature multi-critical point (p*,g*) 
at the classical percolation density p — p* and inter- layer coupling g* ~ 0.16 is demonstrated. The 
quantum critical exponents of the percolating cluster are determined using finite-size scaling. It is 
argued that the associated finite-temperature quantum critical regime extends to zero inter-layer 
coupling and could be relevant for antiferromagnetic cuprates doped with non-magnetic impurities. 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Mg, 75.40.Cx 



Randomly diluted quantum spin systems combine as- 
pects of the percolation problem JlJ with the physics of 
thermal and quantum fluctuations. In systems that can 
be tuned through a T = phase transition as a func- 
tion of some parameter one can hence study divergent 
quantum fluctuations coexisting with classical fluctua- 
tions due to percolation. A multi-critical point, where 
the two types of fluctuations diverge simultaneously, is 
realized in the transverse Ising model with dimensional- 
ity D > 1 |Q, H|. In models with O(N) symmetry and 
N > 2 such a point was believed not to exist, because 
quantum fluctuations were argued to always destroy the 
long-range order on the percolating cluster ||. Several 
studies of diluted 2D Heisenberg antiferromagnets were 
consistent with this scenario ||, k|. However, recent 
quantum Monte Carlo simulations have shown that long- 
range order in the 2D Heisenberg model persists until the 
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and that the percolating clus- 
This implies that the phase 



transition is a classical percolation transition. It also sug- 
gests that a multi-critical point, at which the percolating 
cluster is quantum critical, could be reached by includ- 
ing other interactions. In this Letter it will be shown 
that the 0(3) multi-critical point can be realized in the 
Heisenberg bilayer with dimer dilution, i.e., where adja- 
cent spins on opposite layers are removed together. This 
system is illustrated in Fig. [j], and a schematic T = 
phase diagram is shown in Fig. |^. In analogy with quan- 
tum critical points in clean 2D Heisenberg antiferromag- 
nets ||Tci| , pd| , one can expect a finite-T universal quantum 
critical scaling regime to extend to couplings well be- 
yond the T — critical coupling <?*, possibly all the way 
to decoupled layers (g — 0). This quantum criticality 
could then be realized in layered antiferromagnets doped 
with non-magnetic impurities. It may already have been 
observed in La2Cui_ x (Zn,Mg) x 04, for which recent neu- 
tron scattering experiments [Q show a correlation length 
divergence roughly consistent with the dynamic exponent 
z « 1.3 extracted here. 

The clean 5=1/2 bilayer Heisenberg model has been 
extensively studied in the past |lj| It undergoes a 
quantum phase transition between an antiferromagnetic 







FIG. 2: Schematic T — phase diagram for the Heisenberg 
bilayer with coupling g and a fraction p of the inter-plane 
dimers removed. Percolation and quantum phase transitions 
are indicated by the solid horizontal line and dashed curve, 
respectively. The circle indicates the multi-critical point. 



and a quantum disordered state as a function of the ratio 
g = J2/J1 of the inter- and intra-plane couplings. The 
critical coupling g c w 2.52 and the exponents are consis- 
tent with the expected Jul classical 3D Heisenberg uni- 
versality class. If the system is diluted by randomly re- 
moving single spins, the quantum phase transition is de- 
stroyed because moments are induced around the "holes" 
in the gapped phase. These localized moments order an- 
tiferromagnetically for all g |p"5| . In order to circum- 
vent this "order from disorder" phenomenon, inter-plane 
dimer dilution will be considered here. A dimer is not as- 
sociated with moment formation and hence there is a spin 
gap for large g at any dilution fraction p. When g = 
the system corresponds to two independent site-diluted 
Heisenberg layers. In that case, it was recently shown 
that the sublattice magnetization on large clusters at the 
percolation density scales to a non-zero value in the limit 
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of infinite cluster size ||, |), implying a classical perco- 
lation transition. In the bilayer the long-range order on 
the percolating cluster can then be expected to survive 
up to a critical inter-layer coupling g* > 0, leading to the 
phase diagram shown in Fig. 0. 

In order to extract the multi-critical coupling g* , quan- 
tum Monte Carlo simulations similar to those discussed 
m Ref. | were carried out at the percolation point (p* « 
0.4072538 [|l6)). Two types of boundary conditions were 
used. In periodic L x L systems, dimers were removed 
with probability p* and the largest cluster of connected 
dimers was studied. These clusters have a varying num- 
ber TVi of dimers, with (iVi) ~ AL d , where the frac- 
tal dimension d = 91/48 |J and A w 0.67. Clusters 
were also constructed at fixed size N without boundary 
imposed shape restrictions In this case the corre- 
sponding length-scale is R = N l ' d . The two types of 
clusters will be referred to as L x L and fixed- AT, respec- 
tively. A length R = 0.81 • L will sometimes be used 
for the L x L clusters, so that for a given R the average 
size (N±) ~ R x l d is the same as the size of the fixed- 
N clusters. The calculations were carried out using the 
stochastic series expansion method ]l7| . Temperatures 
sufficiently low to obtain ground state properties were 
used for L x L clusters with L up to 64 and fixed- A" clus- 
ters with A* up to 1024. The results were averaged over 
a large number of dilution realizations; from > 10 3 for 
the largest sizes to > 10 5 for smaller sizes. 

For a cluster with A^i dimers (A^ = 2A^ spins), single- 
plane (a — 1) and full-system (a = 2) staggered structure 
factors and susceptibilities are defined as 

i N ° 

Sa(n,n) = — Pij(S?S*), (1) 

a = l 

Xa (n,7r) = — p ij / dr{S*(T)S](0)), (2) 

i,] = l 

where Py = 1 for sites i, j on the same sublattice and — 1 
for sites on different sublattices. Squared sublattice mag- 
netizations are defined as Q (m 2 ) — (3S a (ir,ir)/N a ), 
where () indicates averaging over dilution realizations. 
The two definitions (m 2 ) and (m?,) should extrapolate 
to the same infinite-size value but the finite-size correc- 
tions can be different. In clean systems, it is known that 
the leading size correction is ~ N^ 1 / 2 [[l9| and this was 
found to be the case also for the diluted single layer 
In Fig. ||, (m\) on L x L clusters is graphed ver- 
sus L~ d / 2 ~ (Ni)^ 1 / 2 for several values of the coupling, 
along with the previous g — results. Quadratic fits 
are also shown. Interestingly, the sub-leading correc- 
tions become smaller, i.e., the behavior becomes more 
linear, as g is increased from 0. For g = 0.1 the data for 
L = 20 — 64 can be fitted to a purely linear form, which 
extrapolates to a non-zero value. In the inset of Fig. |^ it 
is shown that the definition (to|) has more curvature but 
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FIG. 3: Finite-size scaling of the single-plane sublattice mag- 
netization on L x L clusters at the percolation density. In the 
inset, both the single-plane (mi) and full-system (1x12) defi- 
nitions of the sublattice magnetization on L x L lattices at 
g — 0.10 are shown along with 7712 for the fixed-A clusters. 

approaches (m 2 ) for the largest cluster sizes. Results for 
(m 2 ) on fixed- A" clusters are also shown to extrapolate 
to the same infinite-size value, with an over-all smaller 
size-correction. For g > 0.2 the extrapolations give nega- 
tive values, indicating that the sublattice magnetization 
vanishes (the fitted forms can here not be correct for very 
large L, as the asymptotic behavior has to be 1/L d if 
(m|) = 0). Hence, the critical coupling 0.1 < g* < 0.2. 

At (g*,p*), both S a (ir, tt) and Xa(jr, ir) should exhibit 
power-law finite-size scaling. In clean systems the quan- 
tum critical scaling forms are S ~ L 1 ^ 11 and \ ~ L 1+Z ~ n , 
where 77 sa 0.03 is the equal-time spin correlation function 
exponent and z = 1 is the dynamic exponent. Since the 
size fluctuates in the case of the Lx L clusters, the statis- 
tical errors in this case are reduced in the size-normalized 
quantities (S a /N a ) and (xa/N a ), which therefore will be 
studied here. If the exponents are defined according to 

(5 Q (7r,7r)/iV a ) ~ i?^, (3) 
(Xa(TT,n)/N a ) - iT*, (4) 

the dynamic exponent can be obtained from the differ- 
ence; z = 7 X — 75. The best over- all scaling behav- 
ior is seen in (S2/N2) for the fixed-A^ clusters. Based 
on this quantity, the multi-critical point is estimated to 
g* = 0.16 ±0.01, and the exponent 7s = -0.90 ±0.01. A 
log-log plot with data for both fixed- A" and L x L clus- 
ters at g = 0.16 is shown in Fig. ^. The L x L data 
have larger corrections to scaling but for the largest clus- 
ters the behavior is completely consistent with the ex- 
ponent extracted from the fixed- A^ data. The power-law 
scaling in (xa/N a ) sets in at larger system sizes, and 
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FIG. 4: Finite-size scaling of the staggered structure factors 
and susceptibilities. The solid lines are fits to results for the 
fixed-iV clusters; the dashed lines have the same slopes but 
are shifted to match the L x L data for large R. 



also in this case the behaviors seen for L x L and fixed- 
N clusters are consistent with each other. The fixed-iV 
clusters again show a wider range of good scaling and 
give 7 X = 0.38 ± 0.03. The dynamic exponent is hence 
z = 1.28 ± 0.02, where correlations between 7s and 7 X 
have been taken into account in the error estimate. 

According to hyper-scaling theory, the temperature 
dependence of the uniform magnetic susceptibility Xu 
should be governed by the dynamic exponent B0| : 



{Xu) ~ T\N 2 




(5) 



Consistency with the z extracted above from T = quan- 
tities can hence be tested in a non-trivial way. Finite- 
temperature calculations were carried out using suffi- 
ciently large Lx L clusters to completely eliminate finite- 
size effects down to T/J\ = 1/256. Results at g = 0.16 
are shown on a log-log plot in Fig. ||. There is indeed 
a significant linear low-temperature regime, where the 
slope is 0.470 ± 0.005. Using d = 91/48 for D in Eq. (§) 
then gives the dynamic exponent z = 1.29 ± 0.01, in full 
agreement with the T — result. 

Another important quantity is the spin stiffness p s . 
It can be calculated in the simulations using the winding 
number fluctuations pl| in systems with periodic bound- 
ary conditions (i.e., using L x L clusters). Hyper-scaling 
predicts p s ~ L 2 ~ D ~ Z §0 



at a T = critical point. 
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FIG. 5: Temperature dependence of the susceptibility at g = 
0.16. Results for different system sizes are shown along with 
a linear fit to the L = 256 data. The small deviations at the 
lowest temperatures indicate that g* is marginally below 0.16. 



At the percolation point one can expect that D should 
not be replaced by the fractal dimension d of the percolat- 
ing cluster, but by its backbone dimensionality db, which 
is significantly smaller (db ~ 1.643 p3|). This is because 
the spin currents wrapping around the periodic clusters, 
in terms of which p s is evaluated, only flow through the 
backbone. Furthermore, the above scaling form applies 
to systems in which the stiffness takes a finite constant 
value in the ordered phase. However, at p* in the single 
layer (and hence for all g < g*) it scales as (p s ) ~ L - */^, 
where t is the conductivity exponent of percolation and 
v = 4/3 is the percolation correlation length exponent 
[|9], The ratio tjv « 0.983 according to recent sim- 
ulations [p3| . In order to account for the "geometric" 
reduction, the following scaling form is tested here: 



(P.) 



-2-d b 



-t/u 



(6) 



Using the value extracted for z above, the exponent 
2 — db — z — t/u « —1.92. Fig. ^ shows data at g = 0.16 
along with this power-law scaling. There are clearly devi- 
ations, but it is also apparent that the numerical results 
are not yet in the asymptotic scaling regime. Signifi- 
cant corrections to the scaling law L~ l l v were also seen 
in the single-layer case 0. A definite test of the con- 
jectured form (^) hence requires calculations for larger 
system sizes. 

To summarize the results, it has been shown that a 
multi-critical point at the percolation threshold is real- 
ized in the dimer-diluted Heisenberg bilayer at a crit- 
ical inter-plane coupling g* = 0.16 ± 0.01. The dy- 
namic exponent z was determined using both T = and 
T > quantities, with both calculations consistent with 
z = 1.28 ± 0.02. Vajk and Greven have recently stud- 
ied the same model at T > ^5|, with results in good 
agreement with those presented here. 
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FIG. 6: Finite-size scaling of the T = spin stiffness at g 
0.16. The line shows the proposed asymptotic behavior. 



In analogy with finite-temperature quantum critical- 
ity in clean 2D antiferromagnets |[o[ [Tl| , there should 
be a significant universal quantum critical regime con- 
trolled by the point (p*,g*). A very interesting ques- 
tion is then whether this universality could be observed 
even in a single diluted layer. The correlation length 
should diverge as T~ 1//z in the quantum critical regime 
Jio| . O . It is intriguing that this form with z ~ 1.4 has 
recently been observed at high temperatures, in simu- 
lations as well as in neutron scattering experiments on 
La 2 Cui_p(Zn,Mg) p 04 Q. However, since the experi- 
mental system at the percolation point corresponds to 
g < g* there should, in the absence of 3D effects, be a 
T — > cross-over to a different dynamic exponent char- 
acterizing the line (p = p* , g < g*). Although the per- 
olating cluster is ordered, the exponential "renormalized 
classical" behavior [flo[ ^6) cannot apply here because 

Pe = 0. 

The T = dynamic exponent for (p = p*,g < g*) 
is governed by the properties of the percolating cluster. 
Since it is long-range ordered the structure factor expo- 
nent 75 = in Eq. (^). In a clean D-dimensional sys- 
tem with long-range order the staggered susceptibility 
diverges as L 2D K|, and if this holds also at the perco- 
lation point it would imply 7 X = d in Eq. (^) and hence 
z(g < g*) — d = 91/48. This is in contrast to the trans- 
verse Ising model, where there is activated scaling, i.e., 
z = oo at the percolation transition || . The bilayer sim- 
ulations are consistent with z(g < 0.1) = d. Closer to g* 
cross-over effects make it hard to verify this asymptotic 
behavior. 

It should be pointed out that the bilayer coupling used 
here to realize a multi-critical point {p*,g*) is only one 
way to achieve this universality class. In Zn and Mg 
doped cuprates there may be interactions driving the sys- 
tem from g = closer to such a point in an extended 
parameter space. The fact that the sublattice magne- 



tization measured experimentally E3] falls significantly 
below the calculated curve |J as p — > p* indicates that 
such couplings indeed are present. 
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